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Abstract 

Recently lots of efforts have been made to obtain the next to leading order and Landau-Pomeranchuk-Migdal 
corrections to the thermal dilepton emission rate in perturbative QCD. Here we apply these results to the plasma 
created in heavy ion collisions and see wether these corrections improve the comparison between theoretical 
calculations and experimental results for the invariant mass dependence of the dilepton emission rate. In 
particular, we simulate the quark-gluon plasma produced at RHIC and LHC using a 2+1-dimensional viscous 
hydro model. We compare our results to STAR experiment and comment on the need for a non-perturbative 
determination of the dilepton rate at low invariant mass. 


1 Introduction 

The theoretical study of the quark-gluon plasma (QGP) 
[1, 2] is supported by experiments based on ultra- 
relativistic heavy ion collisions (HIC): gold ions are 
used at RHIC (Relativistic Heavy Ion Collider) at the 
Brookhaven National Laboratory and lead ions are 
used at LHC (Large Hadron Collider) at CERN. 

Excellent probes to study the QCP are photons and 
dileptons pairs [3], in fact, as they interact electro- 
magnetically, they have a small cross section w^ith the 
strongly interacting matter inside the QGP. Thus they 
leave the QGP and reach the detector, carrying infor¬ 
mation from deep into the plasma [4, 5]. Moreover w^e 
prefer dilepton to photons for two reasons: photons 
are produced from a bigger background of decays [6] 
while, on the other hand, leptons have a non-null in¬ 
variant mass M, which helps in disentangling various 
sources [7]. 

In fact, the dilepton background is also not small: 
dileptons are produced in every phase of the HIC and 
in several types of processes [7]. Here we are interested 
in thermal dileptons [7] produced by the partonic in¬ 
teractions during the hydrodynamical expansion; these 
dileptons can tell us about the QGP properties . Ther¬ 
mal dilepton are produced mainly in quark-antiquark 
annihilation and Compton scattering and their contri¬ 
bution to the spectrum is important in the intermediate 
invariant mass range, M e [0.2,2.5] GeV. 

The first type of background we encounter consists 
in hadronic reactions at early times. They consist in 
jet-dilepton conversion from the initial hadronic scat¬ 
tering processes and from the photoproduction pro¬ 
cesses. These are hard processes that contribute to the 
spectrum in the high mass range (M>3 GeV). 


Secondly particle decays imprint broad peaks in the 
spectrum, for instance, the contribution from the decay 
of open charm cc e^e~X is also very important in the 
intermediate mass range [9, 10]. In the low mass range 
0.6 < M < 1.1 GeV the decays of vector mesons, i.e. p , 
CO, (p, give a sizeable contribution to the invariant mass 
spectrum [5, 9, 10]. Einally below M < 0.2 pions decays 
from the hadronic phase dominate. 

A big effort has been done to study the thermal 
dilepton production from the QGP in perturbative 
QCD: in particular references [3, 11] provide a complete 
discussion of the leading order (LO) with inclusion of 
anisotropic corrections and [6, 12] supply the passage 
from the LO to the NLO and Landau-Pomeranchuk- 
Migdal corrections. 

In this work we investigate whether higher or¬ 
der corrections in perturbation theory can improve 
the comparison between theoretical predictions for the 
thermal dilepton emission and experimental results. 

In section 2, we introduce the theoretical back¬ 
ground, in particular how to compute the LO dilep¬ 
ton emission rate per unit of 4-volume and per unit of 
4-momentum, what is the effect of NLO and LPM cor¬ 
rections to it and how we compute the invariant mass 
spectrum. In section 3, we explain how we describe 
the hydrodynamical plasma evolution using SONIC 
and the details of the numerical computation for the 
invariant mass spectrum. In section 4, we show our 
results and compare them to experimental result from 
the STAR experiment at RHIC [5]. We then make anal¬ 
ogous computations for LHC and conclude in section 
5. 
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2 Dileptons in heavy ion collisions 

We recall here the perturbative QCD calculations of the 
dilepton rate and explain how to use them in the geom¬ 
etry of heavy ion collisions. We use natural units when 
not stated otherwise and metric signature is In 

perturbation theory, two distinct expansions are made, 
one in the electromagnetic coupling, where the LO is 
sufficient and a second one in the strong coupling. In 
this work we discuss the validity of this second ex¬ 
pansion in the case of the plasma created in heavy ion 
collisions. 

2.1 Leading order dilepton production rate 

The relation between the dilepton emission rate and 
thermal expectation values of electromagnetic current 
correlation function 

Wf'’' = J d*xe-‘‘i'^(P‘(x)r'(0)}, (1) 

is well described in refs. [3, 13, 14]. Here we briefly 
summarize the main results for the specific case of a 
qq e'^e~ process, shown in figure 1. 


Q e- 



the case where the lepton mass is negligible compared 
to the invariant mass M = » m, which is a good 

approximation for electrons, we get [12]: 


d^P 



T 1 

2n^ p - 1 


(4) 


where E = po and p = p^. 


2.2 NLO corrections to the spectra 

As emphasised before, large corrections to the leading 
order dilepton rate arises. The next-to-leading order 
(NLO) is suppressed only by but diverges in the 
small invariant mass limit M ^ 0 (some representative 
diagrams for the NLO are shown in fig. 2). In this soft 
limit the correct result can only be obtained by per¬ 
forming the Landau-Pomeranchuk-Migdal (LPM) re¬ 
summation. The LPM resummation takes into account 
the destructive interference effect between the prompt 
emitted photons and it is summarised by the Feynman 
diagram in figure 3. 


e 




Figure 1: LO Feynmann diagram for the thermal dilep¬ 
ton production from the QGP 


Figure 2: Examples of NLO Feynman diagram for the 
thermal dilepton production from the QGP. 


The number of dileptons produced per unit volume 
and emitted at a given total momentum P = (p^.p^) can 
be expressed trough the dilepton rate R: 

dN^^(x,P) ^ dR(x,P) 
d^xd^P d'^P 



which in turn can be calculated form the quark current 
correlator lVf'’'(P) as: 


dR^'^ 

'd^ 



247l3p2 



xe(P^-4:m^)w^,{P). 


-w) 


(3) 


where m is the mass of the emitted leptons and Q/, i = 
the charges of the Uf massless quarks present 
in the plasma. The strong coupling only enters in the 
quark current correlator W^^(P), which is calculated 
below at leading order but which receives large higher 
order corrections. If we restrict to leading order and to 


Figure 3: Feynman diagram that summarizes the LPM 
effect: in the very dense QGP a quark is rescattered 
many times before it annihilates with its antiparticle. 

In this section we investigate how the NLO and 
LPM corrections contribute to the dilepton spectra. 
To show the results of NLO and higher order correc¬ 
tions, we used the data furnished by Ghisoiu and Laine 
[12] (NLO and LPM to LO) and Ghiglieri and Moore 
(for the LPM at NLO) [6]. It consists in a database for 
the electron-positron and the muon-antimuon emission 
rate as a function of invariant mass M, temperature T 
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and modulus of 3-momentum P for the NLO and for longitudinal flow: 
the LPM corrections (available at [15]). 

= (t cosh C, Xj_, T sinh C), (6) 



M/T 

Figure 4: dielectron emission rate computed at lead¬ 
ing order, NLO [12] , and NLO plus LPM corrections 
[6] as a function of the dimension-less quantity M/T 
for the values P = 0.53 GeV, T = 0.18 GeV and the 
quarks u, d, s have been considered. These corrections 
dominate at small total momentum P. 


so that one can rewrite the differential measure of 
space-time as d^X = idi dC d^Xj _. 

The longitudinal Bjorken flow [18] velocity is de¬ 
fined simply as the distance covered over proper time: 

u^i{T,y) = 7B(l,0,0,t>z) = 7 b(1,0,0, ^), (7) 

where vn = , ^ On the other hand, we obtain 

' Vi-h/T)^ 

the temperature in the C = 0 slice, using a (2+1)- 
dimensional hydrodynamic simulation (see section 
3.1). 

2.4 Dilepton spectra in HIC 

One can rewrite formula (4) using a new parametri- 
sation of the 4-momentum of the virtual photon [11], 
which is better suited for the geometry of HIC: 

= (m_i_ cosh y, pj_ cos (pp, pj_ sin pp, mj_ sinh y^. (8) 


Figure 4 shows a comparison between the dielec¬ 
tron emission rate at LO computed with formula (4) 
and the same with NLO and NLO plus LPM correc¬ 
tions. 

We can notice that the NLO and LPM contribution 
to the emission rate are important for small values of 
the ratio of the invariant mass over temperature M/T. 
Even if the LPM corrections suppress the rate at very 
small M, they increase the rate for intermediate values 
of the invariant mass. 


In the formula above, pp denote the azimutal angle, the 

transverse mass is mj_ = pT + and y the momen¬ 
tum space rapidity 


y = 


1 

Iri — -. 

2 p^ — p^ 


( 9 ) 


Using (8), we write the differential 4-momentum as 
d^P = MdMdyp±dpj_dpp. 

We are now ready to compute the invariant mass 
and rapidity differential spectra: 


2.3 Geometry of HIC and hydrodynamics 
of the plasma 


MdMdy 


rvT ^ 


Over all this work, we will keep separated the longi¬ 
tudinal dynamics (along the collision axis), from trans¬ 
verse dynamics, described by the (2+l)-dimensional 
hydrodynamic model in section 3.1 [16,17]. 

We use the Bjorken model [18, 19] to describe the 
longitudinal expansion and thus we assume the exis¬ 
tence of a "central plateau" structure in the production 
rate of particle as a function of space-time rapidity 


Rapidity C and proper time t = are therefore a 

convenient reparametrization of z and t to describe the 


As we are interested only in the thermal contribution 
to the dilepton emission spectrum, the integration over 
d^X is performed only on the quark-gluon plasma vol¬ 
ume, i.e. on the regions with temperature bigger then 
the critical temperature T(x_l,t) > Tc = 0.17 GeV. 

It is important to note that, in formula (10), the val¬ 
ues of p_L are defined in the LAB reference frame, but 
the emission rate dR/d^P (in formula (4)) has been com¬ 
puted in local rest frame (LRF) of the plasma [11]. While 
integrating over the plasma volume, we have to boost 
the observed LAB frame momenta to the LRF of the 
plasma, which we can plug in the dilepton rate for¬ 
mula (4) computed previously. 
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2.5 Change of frame 

Let us consider the integral over space-time in (10) 
and write it more explicitly using (6): f d^X = 
f TdTdxj_dC Notice that the emission rate (4) de¬ 
pends on x_L and t only through the temperature 
which is given by the hydrodynamical simula¬ 
tion, while the dependence on C is given by the Bjorken 
model, as anticipated. 

At the end of the previous subsection, we noticed 
that, whenever we fix a LAB frame value of and a 
volume element in space in the integral (10), it is neces¬ 
sary to boost the momentum to the LRF of the volume 
element of which we want to compute the emission 
rate. 

In order to do this, we need the boost A^{u^) that 
parametrizes the change of coordinates from the LAB 
frame to the LRF [11]. It is defined by the relative ve¬ 
locity that combines the Bjorken longitudinal 

velocity with the transverse hydrodynamical expan¬ 
sion. Once we found we are able to compute 

(10) considering ^A(^(i/tof)p^^g^. In what follows, 
we show the derivation of u^^^{x^). 

From the hydrodynamical simulations we obtain 
the 4-velocity of the transverse flow (for C = 0), 


Kydro ~ Thydroi'^f'^: 


hydro 


hydro 


/Ol, 


where ykydw = 

measured with respect to the collision axis. However 
the whole system is moving along the longitudinal axis 
with Bjorken velocity Vz, estimated in (7). Thus we will 
need to use the relativistic composition law to find the 
total velocity respect to the LAB frame. 



Figure 5: Scheme of the combination of longitudinal 
and transverse velocity, and the final boost A^. 


The total velocity of a generic element of volume in¬ 
side the QGP is the relativistic sum of Vz and 
v^-Vz^ and Vy = Vz^ ^ 

Then, applying the relativistic addition rule for per¬ 
pendicular velocities Vtot = vi + one ob- 

tains Um(x) = ytot ,yB^xY'\ f ) where yM = 

YB] ■ yhydro- 


3 Numerical simulations 


3.1 SONIC 


The collisions between two heavy ions are not well un¬ 
derstood at early times, before the system thermalizes. 
As soon as the QGP is formed (at proper time Tin ~ 0.5 
fm), its space time evolution is described by hydrody¬ 
namic models [20]. 

We simulate the hydrodynamic evolution of the 
QGP using the software SONIC (Super hybrid mOdel 
simulatioN for relativistic heavy-Ion Collisions), de¬ 
veloped by Romatschke, Luzum and others (the code 
is available at [21]) [16, 17,22, 23,24, 25,26, 27]. In this 
section we summarize the model on which SONIC is 
based on. 

It consists into a (2+l)-dimensional model that takes 
into account only the slice at rapidity C = 0, in which 
the center of mass lies. It combines the pre-equilibrium 
flow, modelled as in Ref. [16], with the hydrodynamic 
phase and the latter with the final hadronisation (which 
does not concern this work). 

SONIC simulates the highly boosted and Lorentz 
contracted nuclei starting with their energy density, 
Ttt = 6(f+ z)Ta(x^), where the function Ta has the shape 
[17]: 


Ta - ^0 



1 + ^-( 


( 11 ) 


R and a are the charge radius and the skin depth param¬ 
eters (the values of these parameters can be found in 
table 1 in reference [17]) and Cq is a normalization con¬ 
stant that controls the final charged multiplicity and it 
is set to reproduce the available experimental data 
The pre-equilibrium radial flow velocity is esti¬ 
mated numerically in [16]: 


vf{T,xA = -^dilnTlixA, (12) 


^The sum of relativistic velocities is not commutative. According to Bjorken model, we make the hypothesis that, in the central rapidity 
plateau, the transverse evolution of the plasma is the same at any rapidity. Thus a generic slice of plasma at a generic rapidity Co evolves 
radially exactly like the Co = 0 slice described by hydro model. We want to obtain the total velocity of a generic element of volume of the 
plasma respect to the LAB frame where pB is measured. The LAB frame, in our case, corresponds to the CM frame, thus every time we fix pB 
in the LAB frame we first have to make a boost along the longitudinal axis to the center of the "Co"pLsma slice" and then add the transverse 
velocity relative to that point, obtained with the hydrodynamics simulations. 
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where t = - z^. The initial energy profile is set to 

be 

e(T,xJ = T2(xJ. (13) 

SONIC includes the relativistic viscous hydrody¬ 
namics solver, (here we use version 1.7), that imple¬ 
ments the evolution of the system using the energy 
density from equation (13) and the flow profile from 
equation ( 12 ). 

The main parameters that need to be set in the 
hydrodynamical simulation are: the freeze-out tem¬ 
perature Tc =0.17 GeV; the initial central temperature 
T =0.37 GeV for RHIC and T =0.47 GeV for LHC; the 
shear viscosity 77 /s = l/4n. 

3.2 Integration of the dilepton rate 

In this subsection, we enunciate the details of the nu¬ 
merical computation for the dilepton spectra ( 10 ). 

Firstly we must introduce the setting to the SONIC 
simulations. The starting proper time is Tstart = 0-5 fm 
and the temporal lattice spacing is 0.001 fm. The space 
grid (which spans the x-y plane) is made of 139 lattice 
sites (per each dimension), separated hy dx = dy = 1 
GeV“^, it covers the squared area [-13.6,13.6]^ fm^. Ev¬ 
ery 500 time steps, (0.5 fm), it takes a "snapshot”, i.e. it 
writes into data files all the measurables, from which 
we use the temperature and the transverse velocity. 

From the inner to outer integration in (10), we first 
computed the integration over zdzdC with the method 
of parallelepipeds. We integrated C in the half range 
[0,0.9] divided in 20 steps and then double the result 
(the integral is symmetric for positive and negative val¬ 
ues of Q. The other integrals (over the transverse co¬ 
ordinates X and y) are computed separately with the 
method of trapezes on the same lattice than the SONIC 
simulation. 

The integral f d(l)p is computed in [0, n / 2 ] with 4 
steps and then we multiply the result times 4, (the sys¬ 
tem is symmetric under rotation with period n/2). The 
extremes of integration for pj_ for RHIC simulations 
have been chosen as in STAR: pj_ e [0.2,15] GeV (we 
note that contributions from pj_ > 15 are negligible), 
while we choose pj_ e [0,15] GeV for LHC and the 
integral is computed in 34 steps. 

The vales of were tabulated in advance as a 
function of T, M and p_i and the values of M that we 
plot always correspond to nodes of the 3D mesh on 
which is tabulated. 

The same was not possible for pj_ because the boost 
shifts its value. Thus for a give p± and T we find the 
corresponding by bi-linear interpolation [17]. 


4 Results 

4.1 Predictions for RHIC 



Figure 6 : Invariant mass spectra for the dileptons emis¬ 
sion at RHIC computed at NLO with LPM correc¬ 
tions for different values of the impact parameter b, 
for central rapidity y = 0 , and transverse momentum 
px e [0.2,15] GeV. 


In this section, we present the results for the mass di¬ 
electron spectrum of the quark gluon plasma ( 10 ) for 
Au-Au collisions at = 200 GeV RHIC and for Pb-Pb 
collisions at Vs = 2.76 TeV at LHC. Moreover we com¬ 
pare the RHIC results with the experimental data from 
the STAR experiment. 

Figure 6 shows our results for the thermal dilepton 
spectrum (10) as function of invariant mass M, for dif¬ 
ferent values of impact parameter b and for LO, NLO 
and NLO + LPM approximation. 

STAR measures the electron-positron pairs from 
Au-Au ions collisions at V^ = 200 GeV, as a function of 
the invariant mass M of the virtual photon and its trans¬ 
verse momentum pj_ [5]. STAR can capture emitted lep¬ 
tons at all azimuthal angles and for momentum-space 
rapidity values |y| < 1. To make a good comparison 
to experiment, we also integrate the dilepton spectrum 
over momentum-space rapidity |y| < 1 . 

Moreover, the data are classified by centrality 
ranges, thus it was necessary to average our results 
over the impact parameter b in order to reproduce the 
centrality ranges. To this aim, we integrated over dif- 
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ferent impact parameters b as ^ 


—{%centrality) 


dbb^(b) 

Jb.i„ 

C’"“^dbb 

'Jumin 


(14) 


where hmin and hmax can be found as a function of cen¬ 
trality in [28, 29] Figure 7 shows the dilepton spec¬ 
trum for RHIC averaged on all the centralities 0-80%, 
which, for Au-Au collisions, corresponds to & = 0 - 13 
fm, for different orders in perturbation theory. 



Figure 7: Invariant mass spectrum for the thermal emit¬ 
ted dileptons computed at LO, NLO and NLO+LPM 
corrections for the full detected centrality range (0-80%) 
and comparison with the relative data from STAR. Er¬ 
ror bars for the star data can be found in [5]. 


Of course, the experimental data from [5] include, 
in addition to the thermal dileptons, high energetic 
electron-positron pairs from pre-equilibrium processes 
and mainly the ones generated in the following decays: 
cx) e^e~n^, e^e~y, 77 ^ e'^e~y, rf, co e^e~, 

,p ^ e"e~ (p e'^e~ , J/p e^e~X . For large invariant 
mass, the thermal contribution quickly becomes small 
in comparison to the other contributions discussed in 
the introduction [8, 7, 9, 10]. In the high mass range 
2 < M < 3 GeV, the contribution from particle decays is 
much more important then the thermal one. For M > 3 
GeV the main contribution to the spectrum is given by 
dielectrons couples produced in pre-equilibrium Drell- 
Yan processes, and our forecast are, of course, much 
smaller then experimental data. 

The region in which the thermal contribution is 
dominant is indeed very small, i.e. roughly 0.2 < M < 
1.5 GeV, up to the p,(jo and p peaks. In this region the 
agreement can be tested. 



M[GeV] 




M[GeV] 


Figure 8: From up to down: Invariant mass spectra 
for the dileptons emission at RHIC computed at LO, 
NLO and NLO with LPM corrections respectively, for 
different ranges of centrality. In each plot the predic¬ 
tions based on (10) are compared with the relative 
datas from STAR. The predicted spectra consider the 
rapidity range |y| < 1 and transverse momentum range 
px e [0.2,15] GeV. 

Figure 8 shows the comparison between our predic- 


^The integral over y and the one at the numerator of (14) has been computed with Simpson's rule for parabolic integration. Considering 
the dependence of the emission rate on the impact parameter (see references [28, 29]), this should give the exact solution of the integral in 
the range b e [0,13] fm. 

^Analogue tables are in ref. [30] for LHC. 
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tions and the experimental data for different range of 
centrality, computed as in (14). Surprisingly, the LPM 
corrections overestimate the number of emitted ther¬ 
mal dileptons at small invariant mass and the NLO 
approximation is the closest to experimental data. In 
the very low mass range M < 0.5 GeV, perturbation 
theory breaks down, and a different approach should 
be studied, for example lattice simulations. 

For M > 0.7 GeV, our results are a bit smaller 
than experimental data as expected, since we do not 
consider all the contributions that are included to the 
experimental data. Moreover we can notice that, the 
more the impact parameter b increases, the farthest are 
the predictions with the experimental data. This is ex¬ 
pected, as for large b, the volume of the plasma pro¬ 
duced is smaller and so the thermal contribution to the 
invariant mass dilepton spectrum is less important. 


4.2 Predictions for LHC 

Figure 10 is the analogue of figure 7 but for Pb-Pb col¬ 
lisions at = 2.71 TeV at LHC, with the difference 
that we kept the rapidity y = 0 and we integrated over 
the full transverse momentum range. The number of 
emitted dielectrons pairs is of course bigger. 



Figure 9: Invariant mass spectra for the dileptons emis¬ 
sion at LHC computed at NLO with LPM corrections 
for different values of the impact parameter b. The 
curves are computed at mid-rapidity y = 0, for the 
transverse momentum range: p± e [0,15] GeV. 


Figure 9 and figure 11 are the analogues of figure 6.3 
(NLO + LPM corrections) and 8.2 (NLO) respectively. 



Figure 10: Invariant mass spectra for the dileptons 
emission at LHC computed at LO, NLO and NLO with 
LPM corrections for collisions with centrality in range 
0-80%. The curves are computed at mid-rapidity y = 0, 
for the transverse momentum range: pj_ e [0,15] GeV. 


5 Conclusion 


The main result of this work is the comparison be¬ 
tween the thermal dilepton rate calculated at LO, NLO 
and NLO+LPM. The higher order corrections becomes 
important for small value of the invariant mass of the 
virtual photon M. Comparing our results with exper¬ 
imental data from the STAR experiment at RHIC, we 
see that for small invariant mass, the NLO+LPM rate 
seems to overshoot the data. This shows that the LPM 
effect, even if it damps the rate at very small M actually 
enhance the rate for M ~ 0.5GeV too much to be com¬ 
patible with the STAR experiment. The NLO seems to 
fit experiment best but overshoot the data for M < 0.5 
GeV. For such small values of M it seems that a non- 
perturbative determination is a must. 

We also performed predictions for the LHC, where 
the plasma phase might become more important in 
comparison to other sources. Results from LHC would 
be very useful to settle the tension between STAR data 
and higher order calculations of the dilepton rate. 
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Figure 11: From up to down: invariant mass spectra for 
the dileptons emission at LHC computed at LO, NLO 
and NLO +LPM corrections for different ranges of cen¬ 
trality. The curves are computed at mid-rapidity y = 0, 
for the transverse momentum range: p± e [0,15] GeV. 
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